Recording particle beam X-ray data

专利摘要:
An apparatus for aligning a set of volumetric data on a particle beam X-ray is described. The apparatus includes a particle source (101) for applying a particle beam directed at an object, the particle source being configured to scan at least a portion of the object. The apparatus includes a particle detector (102) for measuring full depth dose distributions of the respective filiform beam impacts downstream of the object. The apparatus includes a controller (104) configured to calculate reference depth dose distributions corresponding to the respective threaded beam shots based on the stop power information and an assumed position of the object such as represented by the set of volumetric data relative to the respective filiform beam shots, compare the measured depth dose distributions with their corresponding reference depth dose distributions, and calculate the position of the object with respect to the beam shots. filiform on the basis of a result of the comparison.
公开号:BE1025557B1
申请号:E2017/5977
申请日:2017-12-22
公开日:2019-04-15
发明作者:Stappen François Michel P. Vander；Sylvain Jean G. Deffet
申请人:Ion Beam Applications S.A.；Universite Catholique De Louvain；
IPC主号:

专利说明:

Recording of particle beam radiography data
Field of the invention
The present invention relates to the recording of particle beam radiography data. More particularly, the present invention relates to the recording of particle radiography measurements on a volumetric data set, such as a CT data set or an MRI data set.
Prior art
Charged particles (i.e. protons, ions such as carbon ions) have physical advantages over X-rays or gamma rays in the field of radiotherapy. For example, protons of a given energy (that is to say forming a mono-energetic proton beam), have a certain depth of penetration into a target object and do not penetrate beyond this depth, and in addition, they deposit their maximum quantity of energy or dose in the Bragg peak, which corresponds to said penetration depth, that is to say the point of maximum penetration of the radiation in the target volume. The position of the Bragg peak is also associated with the "beam range", which is usually defined as the position where the dose is 80% of the value at the Bragg peak. As the position of Bragg peak depends on the energy of the particle beam, it is obvious that by regulating and modifying the energy precisely, one can place the Bragg peak at a given depth of a predetermined region, such than a tumor, so as to administer the greatest radiation energy to selected points and to preserve the healthy tissue surrounding said points.
The location of the Bragg peak in the tissue must be precisely known, since the critical tissue located near the target tumor could receive an excessive dose, whereas on the contrary, the target tumor could receive an insufficient dose.
US Patent 8,461,559 describes a method for evaluating radiation model data in particle beam radiation applications, in particular in proton beam therapy of a determined target volume of malignant tissue in a patient, comprising the following steps: a) obtaining diagnostic data for a determined target volume to be irradiated; b) calculating a range
BE2017 / 5977 of particles in the predetermined target volume on the basis of the diagnostic data for the determined target volume; c) designing a radiation model with particle beam characteristics based on the calculated particle range and optionally on a calculated dose depth distribution; d) applying a single threadlike beam shot at the target volume determined at a high beam energy compared to the particle beam characteristics of the radiation model; e) measuring the beam range of the single impact of the filiform beam downstream of the determined target volume; and f) comparing the measured beam range with a reference beam range calculated based on the radiation model.
The article “Patient-specific stopping power calibration for proton therapy planning based on single-detector proton radiography” The article “Patient-specific stopping power calibration for planning proton therapy based on single detector proton radiography ), PJ Doolan, M. Testa, G. Sharp, EH Bentefour, G. Royle and H.-M. Lu, Physics in Medicine and Biology, volume 60, number 5, 2015, describes an optimization device that can produce patient-specific calibration curves to convert x-ray computed tomography (CT) data to powers of 'relative stop (HU-RSP) for planning a proton therapy treatment. The difference between a path length equivalent map to digitally reconstructed radiography water (DRRWEPL) through the CT x-ray dataset and a proton radiography (defined as reality) is minimized by optimization of the HU-RSP calibration curve.
Summary of the invention
The present invention aims to provide improved information which can be used to improve treatment with particle radiation.
According to a first aspect of the present invention, the invention relates to an apparatus for aligning a set of volumetric data on radiography data by particle beam according to the first independent claim, comprising:
a particle source for applying a particle beam directed towards an object, the particle source being configured to scan at least a portion of the object with a plurality of filiform beam shots at sufficient beam energy
BE2017 / 5977 high to allow the filiform beam to pass through the object and beyond it, the plurality of filiform beam shots being positioned differently with respect to the object;
a particle detector for measuring dose distributions of full depth of the respective filiform beam shots downstream of the object; and a control device configured for:
calculating information regarding a three-dimensional stopping power map of the object based on a volumetric data set representing the object, calculating reference depth dose distributions corresponding to the respective filiform beam shots based on from the three-dimensional map of stopping power and an assumed position of the object with respect to the respective filiform beam shots, compare the measured depth dose distributions with their corresponding reference depth dose distributions, and calculate a difference between an actual position of the object and the assumed position of the object based on a result of the comparison.
This allows the device to find the actual position of the object relative to each of the filiform beam shots in the volumetric data set with greater accuracy. Thus, the measured depth dose distributions can be mapped to the data in the volumetric dataset with greater precision. This allows, for example, to control the source of particles to irradiate a target region derived from the volumetric data set with greater precision. It also allows the object to be positioned more precisely in a particle radiation system before any therapeutic treatment with particles.
The controller can be configured to, after the controller has calculated the position of the object, adjust the stopping powers of the three-dimensional stopping power map based on at least one of measured depth dose distributions. In this way, the accuracy of the three-dimensional map can be improved. For example, after the volumetric dataset has been aligned with the filiform beam shots with greater precision using the
BE2017 / 5977 calculated position of the object, the estimate of the stopping power can be further improved by again comparing the measured depth dose distribution with the reference depth dose distribution (better aligned).
The controller can be configured to calculate information about the three-dimensional stop power map based on a calibration curve matching voxel values from the volumetric data set with power values d stop, and adjust the calibration curve on the basis of the calculated position of the object with respect to the filiform beam shots and at least one of the measured depth dose distributions. This allows for a more precise determination of the calibration curve, allowing a more precise range estimate to be made, which is an important parameter for treatment planning.
The controller can be configured to repeat the steps of calculating information about the three-dimensional stopping power map, calculating the reference depth dose distributions, comparing the measured depth dose distributions, and calculating the position of the object, based on the adjusted calibration curve. This can lead to an even more precise determination of the position. In addition, the step of adjusting the calibration curve can be repeated using the new position, to obtain an even more precise calibration curve. The entire process can be repeated until the process has converged.
The controller may be configured to, within the framework of the comparison, calculate a cost value based on the differences between the measured full depth dose distributions and the reference full depth dose distributions; and, in the context of calculating the position of the object, finding the position of the object associated with a reduced cost value by repeatedly calculating the dose depth reference distributions associated with positions supposedly adjusted by l object and calculating the corresponding cost value. The calculation of the cost value allows the controller to adjust the assumed position so that the cost function decreases, for example by trial and error, thereby improving the accuracy of the assumed positions. The cost value can indicate how the measured full depth dose distributions look like the reference full depth dose distributions.
BE2017 / 5977
The control device can be configured to perform the calculation of the cost value on the basis of the extraction of curve characteristics, and thus to operate a reduction in dimensionality from each dose of integral reference or measured depth, to extract the cost value. This effectively determines the cost value using the information available in full depth dose distributions.
The controller can be configured to offset the depth of the reference integrated depth dose distribution from the measured integrated depth dose distribution corresponding to a threadlike beam shot, and to calculate the cost value over the basis of the depth offset. This can lead to an even more precise determination of the position of the object.
The calculation of the reference depth dose distributions by the control device can include the simulation of filiform beam shots through the object using at least one of the following methods: launching the rays by formation a plurality of weighted rays; Gaussian nucleus convolutions; a cone collapse; and Monte Carlo simulations. This can lead to an even more precise determination of the position of the object.
The plurality of filiform beam shots applied by the particle source can be transmitted in directions parallel or diverging from the object. Additionally, the controller can be configured to combine the measured full depth dose distributions of the respective filiform beam shots to create a particle x-ray of the object, and to record the particle x-ray with the volumetric data set. representing the object based on the calculated position. This allows for great accuracy throughout the object. Particle radiography can be displayed as an image of the object, for example. It also allows you to combine information from particle radiography with information from the volumetric dataset with great precision, for example in a combined visualization.
The particle source can be configured to apply a plurality of sets of filiform beam shots, the filiform beam shots of each set having a different direction, relative to the object, of the filiform beam shots of the other sets. In addition, the controller can be configured to combine
BE2017 / 5977 the measured full depth dose distributions of the respective filiform beam shots of each set of filiform beam shots to create a particle x-ray of the object corresponding to each set of filiform beam shots, and to record the radiography by particles corresponding to each set of filiform beam shots with the volumetric data set representing the object based on the calculated position. This makes it easier to correct the object’s rotations around all three axes of rotation.
The control device can be configured to determine an affine transformation representing the calculated position of the object with respect to filiform beam shots. An affine transformation, typically comprising rotations and / or translations, is very suitable to cover most of the frequently occurring motion effects.
The volumetric data set may include, for example, an X-ray computed tomography (CT) or an MRI (magnetic resonance imaging) data set.
The particle detector may include at least one of: a multi-layer ionization chamber (MLIC), a multi-layer Earaday cavity, multi-layer scintillation detectors, or a multi-layer stacked dose measurement system. Alternatively, any other particle detection system which is capable of directly or indirectly reconstructing full depth doses can be provided.
According to another aspect of the present invention, a method is proposed for aligning a set of volumetric data on radiography data by particle beam according to the second independent claim, comprising:
obtaining a set of volumetric data representing an object;
the calculation of information representing a three-dimensional map of the stopping power of the object on the basis of a set of volumetric data;
scanning at least a portion of the object with a plurality of filiform beam shots at a sufficiently high beam energy to allow the filiform beam to pass through and beyond the object, the plurality filiform beam shots being positioned differently from the object;
measuring depth dose distributions of filiform beam impacts
BE2017 / 5977 respective downstream of the object;
the calculation of reference depth dose distributions corresponding to the respective filiform beam shots on the basis of the three-dimensional map of the stopping power and of an assumed position of the object with respect to the respective filiform beam shots, the comparing the measured depth dose distributions to their corresponding reference depth dose distributions, and calculating a difference between an actual position of the object and the assumed position of the object based on a result of the comparison.
The comparison may include calculating a cost value based on differences between the measured full depth dose distributions and the reference full depth dose distributions;
and the calculation may consist in finding the position of the object associated with a reduced cost value by repeatedly calculating the reference depth dose distributions associated with a supposedly adjusted position of the object and by calculating the cost value corresponding.
Those skilled in the art will understand that the features described above can be combined in any manner deemed useful. In addition, modifications and variations described with respect to the apparatus may also be applied to the method, and modifications and variations described with respect to the method may also be applied to the apparatus.
Brief description of the drawings
The present invention will be described in more detail below, with reference to the accompanying drawings. Throughout the drawings, similar elements have been indicated with the same reference numbers.
Figure f illustrates an apparatus for aligning a set of volumetric data with particle beam radiography data.
Figure 2 illustrates a perspective view of an apparatus capable of generating and detecting particle radiation.
Figure 3A shows a human head with an indication of a position of a
BE2017 / 5977 filiform bundle.
FIG. 3B represents an IDD measured for the filiform beam indicated in FIG. 3A.
FIG. 4A represents a phantom human head with a general indication of filiform beams.
Figure 4B shows two graphs showing measured IDDs for the filiform bundles shown in Figure 4A.
Figure 5 illustrates a method of aligning a volumetric data set with particle beam measurement data.
Figure 6 shows a graph of a measured IDD and a graph of a simulated IDD.
Figure 7A-F shows graphs of a cost function with respect to optimization parameters.
Figure 8 shows an example of a calibration curve.
detailed description
According to the present invention, a proton radiography (PR) using a multilayer ionization chamber (MLIC) or any other device for measuring the dose of integral depth (IDD) can be used for positioning the patient. Because IDDs carry information about lateral heterogeneities along the path of a particle beam, it is possible to achieve inframillimetric accuracy even with a point spacing of 5 mm. Although this detailed description describes the example of proton radiography in detail, the techniques can be applied to other types of particle radiation, such as alpha particles (helium ions), neon ions, oxygen ions, carbon ions, or any type of beam of particle radiation.
In the following, several terms and acronyms are defined and explained in the context of the present invention. The definitions refer to the terminology used in the context of medical imaging, radiotherapy, particle beam therapy, and in particular proton therapy. CT stands for computed tomography using X-ray data acquired by an X-ray CT scanner. SPR (stopping power ratio) means the ratio of stopping power to water. HU means
BE2017 / 5977
Hounsfïeld units, which are units of the CT data set, typically an indicator of relative electron density. WET (Water Equivalent Thickness) means a thickness equivalent to water, which is the thickness of water that would provide the same built-in stopping power for protons (therefore the same energy loss) as a certain sample of material. For example, a particular trajectory of a proton through an object has an equivalent thickness of associated water. IDD (integral depth dose) means an integral depth dose, which is a dose profile measured along the beam axis (or depth axis) of a proton beam. For each measurement point, the dose is (physically) integrated in the lateral plane relative to the beam axis, hence the name "integrated depth dose". MLIC (Multi-layer ionization chamber) stands for the multilayer ionization chamber, which is a device for measuring an IDD using a stack of ionization chambers. PR stands for proton radiography. A proton X-ray can be acquired, for example, using an MLIC in combination with a spot scan, which provides 2D lateral resolution.
One of the limitations to fully exploit the ballistic benefits of protons in proton therapy is the uncertainty of the exact range within an object, such as a patient. This uncertainty is mainly due to the lack of knowledge of the elementary composition of the tissues of the object. Indeed, treatment planning is based on a conventional CT, which gives information relating to the electronic density of tissues, to which X-rays are sensitive. But the proton range depends on another physical property, the stopping power, which depends not only on the density but also on the atomic composition of matter. Unfortunately, X-ray detection cannot correctly distinguish different materials of similar densities, while protons are specifically sensitive to stopping power. This uncertainty is taken into account up to 1.8% of the range among the 3.5% added as the distal margin.
By imaging on the basis of a proton beam instead of X-rays, it is possible to provide more relevant information regarding the stopping power of the tissues inside the object. However, this presents certain difficulties. First, unlike X-rays, protons do not follow a straight line through matter, but draw sinuous paths due to Coulomb scattering.
BE2017 / 5977 multiple (MCS). Then, the protons pass through thin detectors without delivering all their energies; the pixel intensity response of a simple 2D imager only provides information about the position of the proton, but not about the materials the proton passes through.
Figure 1 shows a diagram of an apparatus for aligning a set of volumetric data with a set of particle beam radiography data. The apparatus includes a particle source 101, a particle detector, and a space for an object 103 between the particle source 101 and the particle detector 102. The apparatus further includes a controller 804, which is configured to control device operations. The control device 804 can be implemented in the form of one or more computer processors, for example. The apparatus includes an 805 memory, which may include volatile and / or non-volatile memory. Memory 805 can include RAM, ROM, FLASH, and / or magnetic disk storage, for example. In memory 805, a program can be stored which, when loaded by the controller 804, causes the controller to control the particle source 801, the particle detector 802 and other components of the apparatus. The control unit 804 can also be configured, via the program stored in the memory 805, to perform calculations and data transformations in order to record the measurement data generated by the particle detector on another set. volumetric data. The device may include a communication port 808 for communicating with other devices, in order to receive sets of data produced with other radiography methods, for example three-dimensional sets of objects representing objects, for example sets of CT data or MRI data sets. These data sets can also be stored in memory 805. The apparatus can further comprise an input 806, such as a user input device for receiving commands from a user, and an output 807, such as an audio display and / or loudspeaker, for outputting data and signals to the user. The particle source 801 can comprise any generator of a particle beam, as is known in the state of the art. For example, the source of particles 801 may include a source of protons which emits a beam of protons in the direction of space for the object 803. The
BE2017 / 5977 particle detector 802 is a detector for particles which are emitted by the particle source 801. For example, the particle detector comprises a proton detector. The detector 802 is configured to detect the location inside the detector where the particle (s) release their energy and the amount of energy released by the particle (s) reaching the detector.
The particle source 801 is configured to apply a particle beam directed towards the space for the object 803. More particularly, the particle source 101 can be configured to transmit a threadlike particle beam. Such a filiform beam is a narrow beam, which can have a finite diameter. Typically, the radiation intensity inside the filamentary beam is highest at the center of the beam and lower near the periphery of the cross section of the beam. For example, the two-dimensional cross-sectional radiation intensity distribution of the filiform beam may be Gaussian. The particle source 801 can be configured to scan at least a portion of the object with a plurality of filiform beam shots at a beam energy high enough to allow the filiform beam to pass through the object and beyond this one. The plurality of filiform beam shots are positioned differently from the object. The filiform beam shots can be applied in sequence, and the object 803 can be moved between the filiform beam shots, so that the filiform beam cuts different parts of the object 803. For example, the apparatus may include a support (205; illustrated in Figure 2) on which the object can rest, and the object can be moved by moving the support. In a particular implementation example, the source of particles 801 and / or the particle detector 802 can be moved and / or rotated. The particle detector 802 is configured to measure dose distributions of full depth of the respective filiform beam shots downstream of the object. Thus, after passing through the object, the particles are received by the particle detector 802 and release their energy inside the particle detector 802. The particle detector 802 is configured to generate an integrated depth dose distribution. In other words, the particle detector 802 detects where, along the direction of the beam, the energy is released by the particles. The energy in planes perpendicular to the direction of the beam is integrated by the particle detector 802. In this way, an integrated depth dose distribution is generated. The control device
BE2017 / 5977 is configured to receive a volumetric data set, for example by means of the communication port 808. As a variant, the volumetric data set is extracted from the memory 805. As a variant, the device comprises a volumetric imager, such as a CT scanner (not shown in the drawing), which generates the volumetric data set by scanning the object. The controller 804 is configured to calculate information regarding a three-dimensional object stopping power map based on a volumetric data set representing the object. This can be done using a transform map or a calibration curve that maps the voxel values from the volumetric dataset to the stopping power values. The transformation map or calibration curve may include a table that maps the voxel values from the volumetric dataset to the stopping power values. The voxel values can for example be CT values or, in the case where the volumetric data set includes an MRI data set, T1 or T2 values. The stopping power can for example be expressed as a stopping power ratio which defines a ratio of the actual stopping power to the stopping power of the water. The transformation card or calibration curve can be stored in the memory 805 or received via the communication port 808. The manner of creating such a transformation card or calibration curve is known to the person skilled in the art. job. Although the following description describes detailed embodiments using a CT data set with CT values (or Hounsfield field units), the techniques can be applied similarly to other types of three-dimensional data sets with voxel values that are expressed in different units, such as MRI data sets, for example T1 or T2 values. Thus, the characteristics described here involving the Hounsfield field units can be modified to process the voxel values produced by the detection modality, such as T1 or T2 values.
The device is configured to detect and store position information, that is, the relative positions of the different filamentary beams relative to the object. In other words, the distance and rotation values between the different filiform beams are recorded with precision. For example, the support may have an actuator which can move the support carrying the object 803 with great precision, for example
BE2017 / 5977 example using a stepper motor. In addition, if the gantry of the device is rotatable, the angles of rotation of the particle source 801 and the particle detector 802 between the filiform beam shots are recorded. Thus, the main imprecision concerning the positioning of the filiform beam shots is due to the imprecision of the position of the object 803 (and of possible deformations of shape of the object itself). The controller 804 can be configured to receive information regarding the position of the object in space for the object 803, and also information regarding the position of the object in the volumetric data set. From this information, the controller 804 can calculate an estimated intersection of each filiform beam with the volumetric data set. Since the relative displacement between the filiform beams is known with precision, after the position of a filiform beam relative to the volumetric data set has been determined, the remaining filiform beams can be located relative to the data set volumetric on the basis of this relative displacement. The controller 104 is configured to associate the appropriate calculated stopping power values from the three-dimensional stopping power map with each filiform beam based on the determined position of each filiform beam. Using these associated calculated stopping power values, the controller 104 calculates reference depth dose distributions corresponding to the respective filiform beam shots. These reference depth dose distributions are thus calculated from the volumetric data set. In this sense, the reference depth dose distributions are simulated depth dose distributions.
The controller is configured to compare the measured depth dose distributions to their corresponding reference depth dose distributions. This can be done in different ways, for example by calculating a cost value that represents the difference between the measured depth and reference dose distributions. The controller 104 is configured to calculate the actual position of the object relative to the filiform beam shots based on a result of the comparison. For this purpose, the difference between the actual position of the object and the assumed position of the object can be calculated. The controller 104 can execute an optimization algorithm by modifying the values for the positions of the filiform beams with respect to the volumetric data set and by recalculating
BE2017 / 5977 the corresponding integrated depth dose distributions, then compare the measured depth dose distributions again with their corresponding reference depth dose distributions. Depending on whether the cost value increases or decreases, the values for the positions of the filiform beam shots can be further modified in the same direction or in another direction, or the process can be stopped if the comparison has reached an acceptable level.
The controller can be configured to calculate the cost value based on differences between the measured full depth dose distributions and the reference full depth dose distributions. These differences can be measured along the depth axis of the full and summed depth dose distributions. The cost values for all the filiform bundles can be summed to arrive at an overall cost value. The control device 104 can be configured to find the position of the object associated with a reduced cost value by repeated calculation of the dose doses of reference depth associated with a supposedly adjusted position of the object and calculate the cost value corresponding.
The control device 804 can be configured to carry out the calculation of the cost value on the basis of the extraction of curve characteristics, and thus to achieve a reduction in dimensionality from each dose of integral reference or measured depth for extract the cost value. The controller 804 can be configured to extract at least one characteristic from each reference integral depth dose or from a measured integral depth dose, and calculate the cost value based on the extracted characteristics.
The controller 804 can be configured to effect a depth shift of the reference integrated depth dose distribution from the measured integrated depth dose distribution corresponding to a threadlike beam shot. The controller can then calculate the cost value based on the determined depth offset.
The control device can be configured to, in the context of calculating the reference depth dose distributions, simulate the shooting of filiform beams through the object using, for example, one or more of the following methods: launching of the rays by forming a plurality of weighted rays; a Gaussian nucleus convolution; cone crushing, and Monte Carlo simulations.
BE2017 / 5977
The particle source 801 and the particle detector 802 can be configured to create a plurality of filiform beam shots that are transmitted in the same direction relative to the object, by means of spot scanning. This can be done, for example, by moving the object in space for the object 803 between successive shots of the filiform beam. The controller 104 can be configured to combine the measured full depth dose distributions of the respective filiform beam shots to create a particle x-ray of the object. For example, particle radiography can consist of collecting the IDDs of the plurality of filiform beam shots. Alternatively, each IDD can be converted to a pixel value. Together, the pixels form an image of the object. Such a conversion can be carried out by a conversion formula, certain nonlimiting examples thereof being the maximum value of the IDD or the depth value of the peak of the IDD. The controller 104 can be configured to record the particle radiography with the volumetric data set representing the object based on the calculated position.
The particle source 801 can be configured to apply a plurality of sets of filiform beam shots, the filiform beam shots in each set corresponding to a particle X-ray having a different orientation relative to the object. Thus, the filiform beam shots of each set are transmitted by the source of particles 801 in a different direction, relative to the object, of the paths of the filiform beam shots of the other sets. This can be achieved, for example, by moving the object 803 between successive shots of the filiform beam. To generate sets with paths having different directions, the apparatus can be configured to rotate object 803 between successive sets of filiform beam shots. In this case, the device may include a rotating robot arm that holds the object. Alternatively, the gantry (not shown) on which the particle source 801 and the particle detector 802 are mounted can be rotated. The controller 804 can be configured to combine the measured full depth dose distributions of the filiform beam shots of a set of filiform beam shots to create a particle x-ray of the object corresponding to that set of beam shots threadlike. This can be done for each set. The controller 804 can be configured to record the x-ray by
BE2017 / 5977 particles corresponding to each set of filiform beam shots, the volumetric data set representing the object on the basis of the calculated position.
In addition, the controller 804 can be configured to determine an affine transformation representing the calculated position of the object relative to the filiform beam shots. In this case, the parameters to be optimized can include the elements of a matrix representing a linear transformation and a three-dimensional vector representing a translation. An affine transformation also makes it possible to take into account the deformations of the object during recording.
The controller 804 can be configured to, after the controller has calculated the position of the object, adjust the stopping powers of the three-dimensional stopping power card based on at least one measured depth dose distributions. Based on the adjusted stopping powers and / or the determined actual position of the filiform beams with respect to the object, the controller 804 can adjust a schedule of particle radiation therapy based on the position setting. assumed from the object, adjusted stopping powers and / or adjusted calibration curve. This planning of particle radiotherapy can be adjusted in a known manner in the state of the art, by correcting the positioning and the applied energy of the particles in a therapeutic particle beam emitted by the source of particles 801.
The particle detector 802 may include at least one of a multilayer ionization chamber (MLIC), a multilayer Faraday cavity, multilayer scintillation detectors, or a multilayer stacked dose measurement system.
Figure 2 shows an illustration of some of the components of a proton radiography device in a perspective view. The same reference numerals have been used in Figure 1 and Figure 2 to indicate similar components. The figure shows the source of particles 801 and the particle detector 802 on opposite sides of an optional support 205. The object 803 to be measured can rest on the support. The drawing also shows a three-dimensional coordinate system 207, with x, y and z axes. In the drawing, the z axis is generally parallel to an imaginary line connecting the center of the particle source 801 and the center of the particle detector 802. This z axis is also called the depth axis. The x-axis and the y-axis are orthogonal to the z-axis.
BE2017 / 5977
One possible goal of proton radiography is to improve the three-dimensional stopping power map, which is typically generated from a three-dimensional CT data set. This conversion step is an important step in the proton therapy treatment planning workflow. This conversion from CT values (Hounsfiel field units) to the stopping power values can be carried out by the following steps. First, a CT data set is acquired and a proton X-ray is acquired. The expected IDDs from proton radiography are determined by simulation, using the CT data set and a conversion curve that maps Hounsfiel field units to the shutdown power ratio. Thus, the Hounsfi field values from the CT dataset are converted to stopping power values, and the IDD corresponding to a particular proton beam is calculated by simulation. The actual acquired IDDs and the simulated IDDs are then compared. For example, the error with respect to the thickness equivalent to water (WET error) is determined as the difference between the thickness equivalent to water calculated from the DLI acquired and the thickness equivalent to water corresponding to the simulated DLI. By optimizing the conversion of Hounsfïeld field units to stopping power ratio, which leads to a smaller WET error, a more precise conversion table between Hounsfïeld units and stopping power is obtained, so that a more precise three-dimensional stopping power map can be derived from the CT data set.
Range error calculation is very sensitive to alignment errors in PR and CT data sets. For example, a 1 mm misalignment would lead to large range errors, especially on regions with high heterogeneities, such as the outer limits of the object. Therefore, this misalignment is a significant confounding factor: while looking for range errors due to the CT-SPR card conversion curve, in fact range errors due to patient alignment during PR acquisition are measured.
The techniques of the invention can be used, among other things, for at least two different purposes: First, image recording can be performed and used to separate range errors caused by misalignment from range errors caused by an imperfect conversion from CT to SPR. For example, a first alignment is improved, then the conversion error can be determined using
BE2017 / 5977 data correctly aligned. Second, techniques can be used to help align the object in a proton apparatus, using the CT image as a guide.
The present invention provides techniques for accurately registering a reference image between PR and CT data sets. In a particular example, the method is based on minimizing the residual error between PR acquisitions (acquired IDDs) and simulations performed on the CT data set (simulated IDDs). Said error can be calculated in different ways. An example is the sum of the square error between each simulated IDD and each measured IDD. Said simulations can be performed, for example, using the following steps to generate a simulated IDD:
1. Convert the CT dataset to a three-dimensional SPR dataset.
2. Determine the WET of a path through the 3D SPR dataset (obtained by conversion from the CT dataset).
3. Shift a reference IDD (for example to the highest energy) by the calculated WET.
4. Weight by a two-dimensional Gaussian curve corresponding to the shape of the proton point.
The steps mentioned above will be explained in more detail elsewhere in this document.
In a particular example, the recording itself can be done in three main steps, all of which are based on minimizing the overall WET error:
1. A recording with two degrees of freedom (using translations perpendicular to the direction of the proton beam) (this step is optional).
2. A recording with six degrees of freedom (using a translation in all three directions and a rotation around any one of three orthogonal axes).
3. Optimization of the CT to SPR curve, then reiteration in step 2.
It will be understood that variations in the degrees of freedom of the recording are possible at each stage. For example, in step 2, recording at six degrees of
BE2017 / 5977 freedom can be simplified in recording with five degrees of freedom by means of a translation in two directions (for example the x and y axes illustrated in figure 2) and the rotation around any of three orthogonal axes. The steps mentioned above are described in more detail below.
A ray tracing algorithm can be used to calculate simulated IDDs, for example to simulate IDDs that have been measured with MLIC. The simulation takes into account the size of the beam and the fact that lateral heterogeneities along the beam path introduce a mixture of the range. Such a model allows the development of a recording algorithm with sub-millimetric precision.
To perform the registration, a cost function can be defined to describe how much the IDD produced by simulation differs from the IDD provided by the MLIC. Next, the CT (converted to SPR) is moved relative to the paths used to calculate the simulated IDD and an interpolation method is used. A CT move provides a different simulated IDD than before, since the simulated path crosses a different part of the CT / SPR dataset. The effect of a displacement is generally greater when lateral heterogeneities are present along the path of the beam. The fluence of a cross section of a filiform beam can be approximated by a two-dimensional Gaussian distribution. In the present invention, the fluence is the radiant energy received by a surface per unit of surface, integrated over the irradiation time.
In a practical embodiment, the standard deviation of this Gaussian distribution can be about 3 mm or more. Because of this, the full width at half height (EWHM) can be greater than 8 mm so that an IDD not only contains information about the WET at the center of the filiform bundle but also near this center , as illustrated in FIG. 3. In FIG. 3 A, a position of a filiform beam through a human head is indicated at number 301. As illustrated, the filiform beam crosses a limit of the skull. FIG. 3B illustrates the IDD measured for this filiform beam, with the depth on the horizontal axis in mm and on the vertical axis the dose in arbitrary unit that the filiform beam emitted in the MLIC. The two peaks show that two sets of protons meet different tissues on their way through the patient. Indeed, in this example, part of the proton beam has crossed the skull while the other part of the proton beam has
BE2017 / 5977 encountered only a soft tissue. In Figure 4, the possible impact of shifting the object from the beam path is illustrated. FIG. 4A represents a phantom human head on a support and an indication 401 of a location where two filiform beams, with a lateral displacement of 1 mm, were measured. FIG. 4B shows two graphs showing the IDDs for these two filiform beams, with the depth on the horizontal axis in mm and on the vertical axis the dose in the arbitrary unit that the filiform beam emitted in the MLIC. We see that a remarkable difference is visible between these two IDDs despite the small displacement.
For example, a simulation can be performed based on the assumption that the filiform beams have a perfect Gaussian shape with a constant standard deviation along the entire path. It has been shown that simulations using this assumption correspond very well to the measurements, even in the case of lateral heterogeneities.
FIG. 5 illustrates in more detail a method of aligning a set of volumetric data with radiography data by particle beam.
In step 501, a "base IDD" is determined, corresponding to an equivalent known thickness of water (WET), which could be any predetermined thickness, for example zero. For example, the basic IDD can be determined by measuring a proton beam that only passes through air before entering the MLIC, in addition to the components of the proton apparatus that are also used during acquisitions of proton radiographs. This basic IDD corresponding to a beam passing through the air only is designated here IDDbasic. The measurement can be performed during the manufacturing or configuration stage, and stored in a memory of the device. Alternatively, the basic IDD can be measured shortly before the patient is placed in the measurement area.
In step 503, a volumetric data set is obtained representing a particular object, for example using a CT scanner or an MRI scanner. In the following description, it will be considered that a CT scanner is used, but this can be replaced by another type of volumetric scanner.
In step 507, the IDDs for which the reference IDD is calculated in step 506, are measured by positioning the object in the proton apparatus, and causing the proton beam generator to send a proton beam toward the object and the proton detector, such as MLIC in this example. MLIC measures
BE2017 / 5977 amount of energy released by the protons, as well as the location (along the z axis) where the energy is released. However, the MLIC integrates the energy of all the small beams in the x plane, y orthogonal to the z axis, the z axis, or the depth axis, being the axis parallel to the direction of the beam of protons . Thus, the measured IDD is similar to the reference IDD as will be calculated in step 506. It will be understood that the order of the steps can be modified; for example, the acquisition of the IDDs using the MLIC in step 507 is independent of the calculation of the equivalent water thickness (WET) 505 and of the calculation of the reference IDDs 506. In addition, will understand that steps 505, 506 and 507 can be performed for any desired number of beams and small beams. In addition, the directions of the beams need not be the same. In general, a set of IDDs measured for a set of parallel or divergent proton beams is called a proton X-ray. Proton beams from a proton X-ray can be arranged to cover a particular area of interest (e.g. the entire object or a particular part of the object) by scanning the area using spacing predetermined between the parallel individual proton beams. Such multiple proton radiographs can be acquired to view the object from different directions.
In step 502, a calibration curve is determined. This calibration curve can be extracted from a memory of the device, for example. The calibration curve can be determined by performing appropriate measurements during the manufacture of the device. The way of determining such a calibration curve is known in the state of the art. The calibration curve maps values from a three-dimensional imaging modality, such as Hounsfield field units from a CT data set, to stopping power values, which can be expressed as shutdown power reports. Alternatively, values of a particular type of magnetic resonance data can be converted to stopping power values using a table. These conversion tables can be created as follows. First, we determine what type of tissue corresponds to a particular value of a voxel in the volumetric data set. Then, the stopping power of this type of tissue can be determined by experience. The voxel value is then matched with the stopping power of the corresponding tissue. If several types of tissue have the same voxel value (for example the same value of
BE2017 / 5977
Hounsfïeld) but different stopping powers, a choice can be made for a particular stopping power, for example according to the types of tissue which are supposed to be present in the object; this can lead to inaccuracies in the calibration table.
In step 504, the CT data set is converted to stopping power values using the calibration curve of step 502. In some embodiments, the stopping power values can be expressed as stopping power ratios, i.e. the ratio of the stopping power of the particular fabric to the stopping power of water.
In step 505, the stopping power encountered along the beam path is calculated, by integrating the stopping powers along the path. For example, the stopping power is calculated as thickness equivalent to water (WET), i.e. the stopping power integrated along the beam path is equivalent to the power d 'stop of a water column having a particular thickness. This thickness is a unit to express the stopping power.
The WET along the assumed beam path can be calculated as follows, by combining the WET from a plurality of small beams. A small beam is a mathematical entity corresponding to a trajectory which has a "substantially zero thickness". Alternatively, a small beam can represent a single particle (for example a proton). Assuming that all the small beams that make up the beam are parallel to each other and that the direction of the beam is parallel to a given z axis, and that (x, y, z) forms an orthogonal coordinate system, the WET at the position (x, y) in the plane can be calculated using the following formula:
Here, the spacing (z) represents the size of the voxels in the z direction, N _z denotes the number of voxels in the z direction, and z /, · denotes the values of z encountered along the beam.
In step 506, reference IDDs are calculated by beam simulation. This is based on the water equivalent thickness of each small beam, which was calculated in step 504, the basic IDD, determined in step 501, and the density distribution of the proton beam. , which is represented by a two-dimensional Gaussian curve in
BE2017 / 5977 the present example. Using this information, the IDD corresponding to a simulated small beam for a position (x, y) would be as follows:
IDD (x, y, z) = IDD _basic (z + WET (x, y)
That is, a given WET of a beam path results in an IDD shift in the opposite direction of the beam, since for a larger WET, the proton releases its energy at a position inside the MLIC closest to the surface of the MLIC facing the object and the proton source.
The simulated IDD, or reference IDD, designated by the IDDdirect of a simulated beam can be a superimposition of the IDDs of the small beams that form the beam. Using the Gaussian distribution of the intensity of the particle beam, this superposition leads to the following equation:
NOT
i = l in which G (x, y) is a two-dimensional Gaussian distribution representing the thickness and intensity distribution of the beam, N denoting the number of small beams which are calculated to calculate the reference IDD, and where ( xi, yî), for i = 1,2, ..., N represent positions of the small beams inside the cross section of the beam.
In step 508, the measured IDDs are compared with the corresponding reference IDDs. This comparison can involve any type of appropriate measurement. First, the IDDs of each measured beam can be compared with their corresponding simulated reference IDDs. Then, the results of the comparison of the individual IDDs can be combined to arrive at an overall quality of the correspondence between the entire set of proton beams and the CT data set. For example, the sum of Quadratic Terror between a reference IDD and a measured IDD can be used to compare the measured IDD to the corresponding reference IDD. Alternatively, Range Terror in terms of an offset error can apply to a reference IDD to minimize the sum of quadratic Terror between the reference IDD and the measured IDD can be used. Again, the gap
BE2017 / 5977 type of the absolute error between a reference IDD and a measured IDD can be used. Other examples are provided elsewhere in this description.
The overall cost value can be calculated by adding the error values of individual IDDs, for example. This can also be done as a sum of squares. In addition, the IDDs can be weighted differently in this summing step. For example, beams crossing a heterogeneous part of the CT data set may receive a higher weight than beams crossing a relatively homogeneous part of the CT data set.
In step 509, it is determined whether the comparison of the IDDs is satisfactory. For example, it is determined whether the overall cost value satisfies certain constraints, for example by checking whether the overall cost value is less than a predetermined threshold.
If the overall cost value does not yet satisfy the constraints predetermined in step 509, the transformation parameters are adjusted in step 511, according to an optimization algorithm, then steps 505, 506, 508 and 509 are repeated with the adjusted transformation parameters.
In a particular implementation example, each parameter of a rigid record is optimized one after the other. These parameters include, for example, a rotation around the x axis, a rotation around the y axis, a rotation around the z axis, a translation in the x direction, a translation in the y direction and a translation in the z direction. Using such an approach can be advantageous from a computational expense point of view. With such an approach, most steps do not require transformation (such as resampling) of the 3D CT. For example, in some embodiments, in step 511, only one parameter is adjusted at a time, using an iterative optimization scheme. For example, a particular parameter will be adjusted for a predetermined number of iterations, or until a certain threshold on the cost value has been reached. In the following iterations, the following parameter will be optimized.
The parameters are optimized using an iterative optimization scheme, that is: the IDDs are calculated taking into account the modified transformation parameters, and the optimization scheme determines the way in which the transformation parameters will be changed again, based on the result of the comparison between the corresponding calculated reference IDD and the measured IDD. We will understand that after
BE2017 / 5977 modification of the transformation parameters, the reference IDDs are calculated again, but the measured IDD remains identical. No further action needs to be taken. Alternatively, a multivariate approach can be used for optimization, which could lead to even better precision, but could be slower than optimizing the parameters one by one.
If a plurality of proton radiographs are acquired in different directions, the different proton radiographs acquired in different directions can then be used to optimize different transformation parameters. A proton radiograph is a collection of IDDs corresponding to a plurality of substantially parallel proton beams. For example, if two orthogonal directions are measured, one with proton beams parallel to a first axis and another with proton beams parallel to a second axis, the IDDs of proton radiography with parallel proton beams to the first axis can then be used to optimize the rotation around the first axis, and to optimize the translation along two axes perpendicular to the first axis. The IDDs of proton radiography with proton beams parallel to the second axis can be used to optimize rotation around the second axis, and translation along two axes perpendicular to the second axis. This works even if the acquisition directions are not perpendicular.
To estimate the roll, i.e. rotation about an axis that is orthogonal to both the first axis and the second axis, data from the two proton radiographs may preferably be used. In principle, it is possible to use all the x-ray data to optimize each transformation parameter, however, calculation efficiency is obtained using the selections described above from proton x-rays. This principle can be extended using more proton radiographs in more directions. In this case, more axes of rotation and translation can be used, which can then be reduced to three parameters of rotation and three parameters of translation.
It is observed that, although in the above description, rigid recording parameters are described, the method also extends to more general transformations such as affine transformations or non-rigid general recording methods. The parameters of such transformations can be optimized a
BE2017 / 5977 by one or using multivariate techniques.
After adjusting the transformation parameters in step 511, the method comprises steps 505, 506, 508 and 509, until the comparison value is found sufficient in step 509, after which the method proceeds to l 'step 512.
In an optional step 512, it is determined whether to update the calibration curve using the calculated transformation parameters. If so, the process goes to step 502 to determine the calibration curve using the calculated transformation parameters. For example, after determining the precise position of the CT scan in relation to the filiform beams, precise mapping can be used to recalculate the calibration curve using the measured depth dose distributions and their supposedly adjusted positions. To calculate or improve the calibration curve or the CT to SPR transformation based on the measured IDDs, a method can be used as follows.
The patient's 3D shutdown power map can be adjusted based on proton radiography data by producing a patient-specific calibration curve. Once this curve is obtained, it can be used, for example, to convert the X-ray CT planning into a fitted 3D SPR map.
Such a calibration curve may include a piecewise linear function which maps Hounsfiel field units (HU) to the stopping power against water (SPR). In this approach, the HU and SPR of different materials can be determined. For example, the following table shows the HU and SPR for some well-known tissue types.
Material name HU SPR Air -1024 0 Defibrated lung -741 0.26 Fat 3 -98 0.94 Fat tissue f -55 0.98 Soft tissue 0 1 Skeletal muscle 1 5 1 Skeletal muscle 2 40 1.03
BE2017 / 5977
Skin 1 72 1.05 Connective tissue 100 1.07 Sternum 385 1.18 Humerus 538 1.25 Femur 688 1.35 Skull 999 1.53 Cortical bone 1524 1.82 Titanium 3200 3.19
Board
Table: HU versus SPR for different materials (see Woodard, H. Q., White, D. R., 1986, “The composition of body tissues”, Br. J. Radiol. 59: 1209-18)
Table 1 shows a discrete relationship between HU and SPR for different materials. A continuous relationship can be obtained by assuming that the SPR of these materials is connected by linear segments (i.e. using linear interpolation). This results in a calibration curve such as that shown in Figure 8.
Using the measured DLI information, the calibration curve can be improved. This can be done in several ways. For example, an improved estimate of the SPR of different materials can be calculated. From this information, the improved calibration curve can be determined using linear interpolation, for example. An improved estimate of the SPR of different materials can be calculated, for example, using the following equation:
x = argmin (| 71x - b | ² ), subject to x> 0,
X x being a vector of dimension Μ, M being the number of materials to be evaluated. The components of the vector x designate the SPR associated with each of the M distinct materials. As the Hounsfield field unit of these M materials is known, they can be used to update the above table, and can be used to improve the curve shown in Figure 8. The materials to be used for this purpose can be chosen arbitrarily. For example, the materials represented in the vector x can be chosen so as to correspond to typical human tissues, such as those listed in
BE2017 / 5977 table 1.
In addition, A is an N x M matrix which represents an estimate of the amount of interaction of N particle beam shots with each of the M materials. This amount of interaction can be an indication of a length traversed in a material by a threadlike beam shot. This can be estimated in different ways, for example, based on the CT dataset and the trajectory of each threaded beam shot through the CT dataset (while applying for example the weighting function Gaussian to take into account the intensity of the beam). For each voxel in the CT data set, the material can be determined as material with the HU closest to the HU of the voxel. Alternatively, the material of the voxel is determined as a mixture, for example, of two different materials, namely the material having the smallest HU which is larger than the HU of the voxel and the HU of the material having the largest HU which is smaller than the HU of voxel. For each voxel and for each material, the fraction of the voxel that is occupied by this material can be estimated. Then, the length ay traversed in a material j by a filiform beam i can be calculated as
Oyj = ((estimated length of the filiform beam i through voxel v) all voxels v x (estimated fraction of voxel v which is occupied by material j))
In addition, b is an N-dimensional vector representing the WET values derived from N IDD measured from the N shots of filiform beam. For example, a table can be used that maps the R-80 range of an IDD to the corresponding WET value, taking into account the initial energy of the particles.
With reference to the equation x = argminfllfex - b | ² ), subject to x> 0, additional constraints on the variation of the stopping power of the materials can be defined. For example, equality constraints can be used for certain known materials, such as materials corresponding to air and water. In addition, the interpolation method for creating the calibration curve can be replaced by a variant. For example, linear interpolation can be replaced by another type of interpolation.
By adding more equality constraints, i.e. by imposing that the SPR of certain materials are at predefined values, it is possible to optimize the SPR
BE2017 / 5977 from just a few selected materials, or even just a single material. In particular, this can be useful when a specific material is known to have a considerable impact on the range uncertainty. A typical example is a metal implant. In such a case, for example, a single filamentary beam can be used to optimize the SPR corresponding to this specific material, although several shots can increase the quality of the result.
As on the one hand, the quality of the optimized calibration curve may be impacted by installation errors and, on the other hand, the accuracy of the recording algorithm could be impacted by an incorrect calibration curve , the accuracy of the two methods can be improved by combining them in an iterative process. For example, recording of proton radiography data followed by optimization of the calibration curve can be done iteratively until no significant change in the outputs of the two processes is observed.
Once the calibration curve is complete, the SPR for each voxel can be determined by applying the calibration curve to the Hounsfield field unit of the voxel. In this way, an improved three-dimensional map of stopping power is obtained.
After the calibration curve has been updated in step 502, the CTSPR conversion can be performed again in step 504 using the updated calibration curve. Then, either the process can be terminated, and other operations can be carried out as described for step 510, using the updated calibration curve and / or an updated three-dimensional SPR card. day. As a variant, the process takes place by re-optimizing the positions, using the updated three-dimensional SPR card, in steps 505, 506, 508, 509 and 511.
If the cost value satisfies the constraints predetermined in step 509, and no updating of the calibration curve is decided in step 512, the method proceeds to step 512, in which general operations are performed using the determined position of the object represented by the CT scan with respect to the proton beams (which is relative to the proton apparatus), and / or using the d curve updated calibration or updated three-dimensional SPR card. For example, the information can be used to calculate beam ranges to target a certain region of the object. Another example is a mixed visualization of the whole
BE2017 / 5977 of CT data with proton radiography. Other uses are also possible.
The particle source can be configured to emit a threadlike beam shot having a finite diameter, and therefore comprising a plurality of protons emitted in a narrow beam. The depth dose distribution calculation of steps 505 and 506 described above takes this into account when calculating the reference IDD corresponding to a filiform beam as the weighted sum of a plurality of IDDs caused by individual protons or small bundles. The weights of the weighted sum may represent a beam intensity throughout the cross section of the filiform beam. This approach allows, among other things, to account for the fact that, due to the finite size of the cross section of the filiform beam, multiple Bragg peaks can be detected from a single shot of filiform beam.
Comparing the measured IDD with the reference IDD, represented in step 508, may involve carrying out an extraction of curve characteristic, which achieves a dimensionality reduction of each dose of reference or measured integral depth by extracting a single value or several values which respectively describe the measured IDD and the reference IDD. These extracted characteristics can then be compared. For example, such feature extractors can be based on depth values, such as peak, R90, R80, center of mass, integral below the curve, median, mean, variance, flatness, flattening, or any statistical moment. . Such characteristics can be determined for both the measured IDD and the reference IDD, and then compared. This leads to a comparison of the measured DLI and the reference DLI.
Comparing the measured IDD with the corresponding reference IDD, shown in step 508, may also involve performing a depth shift of the measured IDD or the reference IDD relative to to the other. The depth offset can be determined for the one for which the comparison shows the greatest similarity, for example, the depth offset corresponding to the lowest cost value found. As in other examples, the comparison may involve a sum of differences, a sum of differences squared, variance, among others. The depth offset necessary to obtain the greatest similarity between the curves can also, by itself, be used as a value
BE2017 / 5977 of cost indicative of the similarity of the measured DLI and the reference DLI.
The comparison illustrated in step 508 can also be preceded by a normalization step: first the reference IDD and the measured IDD can be normalized, so that the comparison provides more relevant results. The actual comparison, for example the sum of the differences squared or another, can be performed on the normalized DLIs.
The coordinate systems used for device calculations and settings can be based on the definitions of IEC 61217, for example. The symbols Rx, R _y and R _z are used in the present invention to designate the matrices with three rotations whose axes of rotation are respectively x, y and z. The associated rotation angles are named θ _χ , d _y and θ _ζ . Translations along the x-axis, the y-axis, and the z-axis are designated by the symbols T _x , T _{y and} T _z , respectively. The transformation can be applied according to the following convention:
y ' - R _y R _z R _x xthere + H ¹ ________________________________________1 z 'zT. ^z .
In other embodiments, the transformations can be modified, for example, the order in which the transformations are applied can be modified or the coordinate system can be modified, leading to different values of the transformation parameters. The transformation parameters can be summarized in a correction vector V, as follows:
V = [θχθγθ _Ζ Τ _χ ΤγΤ _Ζ ]
In order to take full advantage of all the information provided by proton radiography with an MLIC, the cost function can take into account range differences at the point position, as well as the range mixture due to the presence of lateral heterogeneities along the beam path. To this end, the use of an IDD simulation process can be integrated into the cost function evaluation. Such a simulation process, capable of accurately modeling the mixture of the range, makes it possible to carry out a high-quality comparison with the measured IDD.
In the following, several examples of suitable cost functions are described, which can be used to determine the cost value used in the process
BE2017 / 5977 optimization.
In the first example, the squared error between the measured IDDs called IDD _mea s and a simulated IDD, called IDD or IDDdirecte, par
SQE (x, y) = f ^l 0 ^detector (lDDdirect (x, y, z) - IDD _meas (x, y, ztf dz, in which Idetector designates the length of the detector in the z direction.
By summing this error over the available measured IDDs and dividing by the number of points (i.e. the number of measured IDDs), a cost function is obtained, named CF, which represents the average error at minimize:
_ 2x ^, y ^ SQE (Xi, yi) ⁿ idd 'in which Nidd denotes the number of IDD measured and the summation Σ _Χ ί, _γ ί SQE (xi, yî) denotes a summation of the SQE for all pairs ( xi, yi), corresponding to the measurement points of the measured IDDs, with i = 1.2, ..., Nidd.
Alternatively, instead of the SQE quadratic error described above, the standard deviation of the difference between the simulation and the measurement can be used:
STD (x, y) = std ( _direct IDD (x, y, z) - IDD _meas (x, y, z)) the standard deviation being calculated from the values of the subtraction for the values of z for which the measure is available. Typically, these values of z are in the range of 0 to Idetector By summing STD (x, y) for the available values of (x, y) and dividing by the number of points, a cost function is obtained, named CF, which represents the value to minimize:
„„ Σχ ,, γ ^ ΤΟ / χι, γι)
CF = —— ------, ⁿ idd in which Nidd denotes the number of IDD measured and the summation Σ _Χ ί, _γ ί STD (xi, yî) denotes a summation of STD for all the pairs (xi, yi), with i = 1, 2, ..., Nidd.
The range error map between simulated IDDs and those being measured is very sensitive to misalignments. Therefore, an average range error can be used to define an appropriate cost function. Due to the mixture of range and the resulting complex form of the IDD, there are different approaches possible to determine this range error. The range error is the difference between the expected range (s) based on the reference DLI and the actual range (s) such as
BE2017 / 5977 represented by (for example the peak or peaks of) the measured DLI. This difference, named RD below, can be evaluated for example by calculating the range offset which minimizes the squared error between the simulated curve and the measured curve:
By summing this RD (x, y) error over the available measured IDDs and dividing by the number of points (i.e. the number of measured IDDs), a cost function is obtained, named CF, which represents the average error to minimize:
ⁿ idd in which Nidd denotes the number of IDD measured and the summation Σχί, γι RD (xi, yj denotes a summation of RD for all pairs (xi, yi), with i = 1, 2, ..., Nidd .
Another cost function given by way of example considers a single value for each Bragg curve or IDD: the position of the maximum dose in the IDD:
M / lX _direct (x, y) = arg max IDD _direct (x, y, z),
Z and
MAX _meas (x, y) = argmaxIDD meas (x, y, z), Z
Different cost functions can be defined on the basis of MAXdirectfx, y) and MAXmeas (x, y). For example, a cost function can be based on the difference squared of these values:
y) MAX _meas (x> y))
AiIDD in which Nidd denotes the number of IDDs measured and the summation Σχ ,, yi denotes a summation for all pairs (xi, yi), with i = 1, 2, ..., Nidd.
Another example of a cost function is based on the mean squared error after correction for the range error, as for example in the following definition of SQE2:
This would lead to the following cost function:
ⁿ idd
BE2017 / 5977 in which Nidd designates the number of IDDs measured and the summation Σ _χ ; _γ ί SQE2 (xi, yj denotes a summation of SQE2 for all pairs (xi, yi), with i = 1, 2, ..., Nidd.
Instead of the mean squared error, we can also consider the variance of the difference between the reference IDD and the IDD measured after correction for the range error, as for example in the following definition of VAR:
VAR (x, y) = Var ( _direct lDD (x, y, z + RD (x, y)) - lDD _meas (x, y, z) j the variance being calculated from the values of the subtraction for the values of z for which measurement data is available. Typically, these values of z are in the range of 0 to Idetector. Instead of the variance, the standard deviation can be used: This would lead to the following cost function _ Hx ^, y ^ VAR (xi, yjj ⁿ idd 'in which Nidd denotes the number of IDD measured and the summation ^ _x ; _y iVAR (xi, yi) denotes a summation of VAR for all pairs (xi, yi), with i = 1, 2, ..., Nidd.
One or more types of optimizations can be used, particularly in steps 509, 511. For example, a method known as "coordinated descent", or a variant thereof, can be used. Such a process benefits from a reduced computation time in the case of small parallel beams and analytical simulations since geometric transforms can be performed directly on the WET 2D (projected) instead of being on the 3D volume. With this method, each of the six parameters (the three angles of rotation and the three translations) is determined sequentially by a minimization algorithm which, in this case, is the golden ratio method, for example. Other optimization methods can also be used, such as a gradient descent method, for example.
Figure 6 shows, by way of example, a graph of a measured IDD of a proton beam through a head phantom, and a graph of a simulated reference IDD 602 based on the same path through the same head ghost. The baseline DLI was calculated, in this case, based on a CT scan of the lead phantom which was converted to a three-dimensional power shutdown reports (SPR) map. The measured DLI was measured using a proton source and a detector including an MLIC. In Figure 6, the horizontal axis represents the "depth" axis in mm, that is, the distance along the ‘z’ axis. The vertical axis is
BE2017 / 5977 the energy detected, in an arbitrary unit. The data measured in this experiment were acquired on a head phantom (CIRS, USA) using Giraffe dosimetry (IBA, Belgium). The beam energy was 210 MeV and the dot spacing was 3 mm. On the other hand, simulated data was obtained according to the equation for
IDDdirect, as described above:
Figure 7 shows an example of how the cost function varies with the optimization parameters θ _χ , d _y , θ _ζ , T _x , T _{y and} T _z . The shape of the cost function with respect to the six parameters (three angles and three translations) determines how a given algorithm converges to the optimal solution. Figure 7 shows the value calculated by the cost function:
Lx j, yjSQE (Xi, yi)
CF = —— ------, ⁿ idd for angles of rotation ranging from -2 to 2 degrees and translations ranging from -5 to 5 mm. Only six graphs are shown, showing different combinations of the six transformation parameters on the two horizontal axes and the value of the cost function CL on the vertical axis. The graphs for other pairs of transformation parameters are assumed to have similar shapes. For example, the shape of the graph of the cost function with respect to the rotation around the x and y axis (Fig. 7A) is similar to the shape of the graph of the cost function with respect to the translation along the axis x and the rotation around the y axis (Fig. 7F). The apparent convex shape of the function allows the use of various well-known optimization tools. Note that the form of the cost function depends on many factors, including the patient's anatomy.
Certain aspects of the invention may be suitable for implementation in the form of software, in particular a computer program product. The computer program product may include a computer program stored on a non-transient computer readable medium. In addition, the computer program can be represented by a signal, such as an optical signal or an electromagnetic signal, carried by a transmission medium such as a fiber optic cable or air. The program
BE2017 / 5977 IT can be partially or entirely in the form of source code, object code, or pseudo code, suitable for execution by a computer system. For example, the code can be executable by one or more control devices or processors.
The examples and embodiments described here serve to illustrate rather than limit the invention. A person skilled in the art can conceive of other embodiments without departing from the spirit and scope of the present invention, as defined. by the appended claims and their equivalents. The reference signs placed in parentheses in the claims should not be interpreted as limiting the scope of the claims. The elements described as separate entities in the claims or description can be implemented as a single hardware or software element combining the features of the elements described.
The present invention has been described above with reference to a number of exemplary embodiments as presented in the drawings. Modifications and variant implementations of certain parts or elements are possible, and are included within the scope of protection as defined in the appended claims.

权利要求:
Claims (12)
[1]
1. Apparatus for aligning a volumetric data set with particle beam radiography data comprising a particle source (101) for applying a particle beam directed towards an object, the particle source being configured to scan at least a part of the object with a plurality of filiform beam shots at a sufficiently high beam energy to allow the filiform beam to pass through and beyond the object, the plurality of filiform beam shots being positioned differently in relation to the object;
a particle detector (102) for measuring dose distributions of full depth of the respective filiform beam shots downstream of the object; and a control device (104) configured to:
calculating information regarding a three-dimensional stopping power map of the object based on a volumetric data set representing the object, calculating reference depth dose distributions corresponding to the respective filiform beam shots based on from the three-dimensional map of stopping power and an assumed position of the object with respect to the respective filiform beam shots, compare the measured depth dose distributions with their corresponding reference depth dose distributions, and calculate a difference between an actual position of the object and the assumed position of the object based on a result of the comparison;
characterized in that the control device (104) is configured to:
as part of the comparison, calculate a cost value based on the differences between the measured full depth dose distributions and the reference full depth dose distributions; and in the context of calculating the position of the object, finding the position of the object associated with a reduced cost value by repeatedly calculating the reference depth dose distributions associated with a supposedly adjusted position of the object and calculating the corresponding cost value;
and in that the control device (104) is configured to extract at least one
BE2017 / 5977 characteristic from each reference integral depth dose or measured integral depth dose, and calculate the cost value on the basis of said at least one extracted characteristic.
[2]
2. Apparatus according to claim 1, in which the control device (804) is configured to, after the control device has calculated the position of the object, adjusting the stopping powers of the three-dimensional power map. stop based on at least one of the depth dose distributions measured.
[3]
The apparatus of claim 1, wherein the controller (804) is configured to calculate information regarding the three-dimensional stop power map based on a calibration curve matching voxel values of the volumetric data set with the stopping power values, and the controller (804) being configured to adjust the calibration curve based on the calculated position of the object relative to the filiform beam shots and at least one of the depth dose distributions measured.
[4]
The apparatus of claim 3, wherein the controller (804) is configured to repeat the steps of calculating information regarding the three-dimensional stopping power map, calculating reference depth dose distributions, comparing the measured depth dose distributions, and calculate the position of the object, based on the adjusted calibration curve.
[5]
5. Apparatus according to claim 1, in which the control device (104) is configured to effect a depth shift of the reference integrated depth dose distribution with respect to the measured integrated depth dose distribution corresponding to a shot. of filiform beam, and to calculate the cost value based on the depth offset.
[6]
6. Apparatus according to claim 1, in which the control device (104) is configured to, in the context of the calculation of the dose doses of reference depth, simulate the shots of filiform beam through the object using at least one
BE2017 / 5977 of the following methods: launching a ray by forming a plurality of weighted rays; a Gaussian nucleus convolution; a cone crush; and Monte Carlo simulations.
[7]
7. Apparatus according to claim 1, the plurality of filiform beam shots defining parallel or divergent paths through the object, and the control device (104) being configured to combine the measured full depth dose distributions of the shots of respective filament beam to create a particle x-ray of the object, and record the particle x-ray with the volumetric data set representing the object based on the calculated position.
[8]
8. Apparatus according to claim 1, the particle source (801) being configured to apply a plurality of sets of filiform beam shots, the filiform beam shots of each set defining paths which have a different direction, relative to the object, paths of the filiform beam shots of the other assemblies, and the control device (804) being configured to:
combine the measured full depth dose distributions of the respective filiform beam shots of each set of filiform beam shots to create a particle x-ray of the object corresponding to each set of filiform beam shots, and record the corresponding particle radiography at each set of filiform beam shots with the volumetric data set representing the object on the basis of the calculated position.
[9]
9. Apparatus according to claim 1, the control device (804) being configured to determine an affine transformation representing the calculated position of the object relative to the filiform beam shots.
[10]
10. Apparatus according to claim 1, the volumetric data set comprising a CT data set or an MRI data set.
BE2017 / 5977
[11]
11. Apparatus according to claim 1, the particle detector (802) comprising at least one of: a multilayer ionization chamber (MLIC), a multilayer Faraday cavity, multilayer scintillation detectors, or a measurement system multilayer stacked dose.
[12]
12. Method for aligning a volumetric data set with particle beam radiography data, comprising:
obtaining a set of volumetric data representing an object (503);
computing information regarding a three-dimensional object stopping power map based on a set of volumetric data (504);
scanning at least a part of the object with a plurality of filiform beam shots at a sufficiently high beam energy to allow the filiform beam to pass through and beyond the object, the plurality of filiform beam shots being positioned differently from the object, and measuring depth dose distributions of the respective filiform beam shots downstream of the object (507);
calculating reference depth dose distributions corresponding to the respective filiform beam shots based on the three-dimensional map of stopping power and an assumed position of the object relative to the respective filiform beam shots (506 );
comparing the measured depth dose distributions to their corresponding reference depth dose distributions (508); and calculating a difference between an actual position of the object and the assumed position of the object based on a result of the comparison (509, 511);
characterized in that the comparison includes calculating a cost value based on the differences between the measured full depth dose distributions and the reference full depth dose distributions; and calculating the position of the object includes finding the position of the object associated with a reduced cost value by repeatedly calculating the reference depth dose distributions associated with a supposedly adjusted position of the object and calculating the corresponding cost value;
BE2017 / 5977 and in that the method comprises extracting at least one characteristic from each reference integral depth dose or measured integral depth dose, and calculating the cost value on the basis of said at least one extracted characteristic.

类似技术:

---

公开号 | 公开日 | 专利标题

---

US9830718B2|2017-11-28|Image processor, image processing method, and treatment system

---

Jia et al.2012|A GPU tool for efficient, accurate, and realistic simulation of cone beam CT projections

---

JP6053792B2|2016-12-27|Integration of user input and modification vector field modification in variable image registration workflow

---

WO2010141583A2|2010-12-09|System and method for dose verification radiotherapy

---

FR2872659A1|2006-01-06|METHOD AND APPARATUS FOR DIRECT RECONSTRUCTION IN TOMOSYNTHESIS IMAGING

---

CN101473348A|2009-07-01|Method and system for error compensation

---

JP5844576B2|2016-01-20|Image generating apparatus, method, and program

---

US8488850B2|2013-07-16|Isotropic resolution image reconstruction

---

Wein et al.2011|Self-calibration of geometric and radiometric parameters for cone-beam computed tomography

---

CN108338802B|2019-07-12|Method for reducing image artifacts

---

EP3160587A1|2017-05-03|Method of producing a radiometric physical phantom of a biological organism and physical phantom produced by this method

---

Cassetta et al.2019|Accuracy of low‐dose proton CT image registration for pretreatment alignment verification in reference to planning proton CT

---

CN109414234A|2019-03-01|System and method for being projected from the 3D data set generation 2D being previously generated

---

EP3145587B1|2018-06-27|Method for estimating the dose administered by an external radiotherapy system

---

BE1025557B1|2019-04-15|Recording particle beam X-ray data

---

Brost et al.2018|A mathematical deconvolution formulation for superficial dose distribution measurement by Cerenkov light dosimetry

---

CN111679311A|2020-09-18|System and method for dose measurement in radiation therapy

---

Koch et al.2008|Automatic coregistration of volumetric images based on implanted fiducial markers

---

EP3153102A1|2017-04-12|Method for determining a dose of radiation applied to a patient

---

Polf et al.2021|Applications of machine learning to improve the clinical viability of Compton camera based in vivo range verification in proton radiotherapy

---

Wu et al.2007|Novel image registration quality evaluator | with an implementation for automated patient positioning in cranial radiation therapy

---

Auer et al.2016|Implementation of a pre-calculated database approach for scatter correction in SPECT

---

Lefkopoulos et al.2001|Aspects physiques et méthodologiques de l’imagerie multimodalités et principes de planification dosimétrique pour la radiothérapie conformationnelle tridimensionnelle

---

Jurkovic2004|Evaluation of a desktop computed radiography system for IMRT dosimetry

---

Ono et al.2015|Resolution improvement of point dose distribution in intensity modulated radiation therapy

同族专利:

---

公开号 | 公开日

---

EP3338860A1|2018-06-27|

---

BE1025557A1|2019-04-08|

引用文献:

---

公开号 | 申请日 | 公开日 | 申请人 | 专利标题

---

---

EP2229981A1|2009-03-17|2010-09-22|Paul Scherrer Institut|A method for evaluating radiation model data in particle beam radiation applications|

---

法律状态:
2019-05-16| FG| Patent granted|Effective date: 20190415 |

---

2020-08-27| MM| Lapsed because of non-payment of the annual fee|Effective date: 20191231 |

优先权:

---

申请号 | 申请日 | 专利标题

---

EP16206519.7|2016-12-22|

---

EP16206519.7A|EP3338860A1|2016-12-22|2016-12-22|Registration of particle beam radiography data|

---